C PROGRAM CCUNBENDK ***********************************************************
C
C               Remember to change version number in first write statement.
C               Original pre-history versions called CCUNBEND and CCUNBENDA. 
C
C       VX 1.0  RH      21.6.84 CCUNBENDB - This went through a variety of
C                               changes before version numbers were started.
C       VX 2.0  JMB     20.3.86 Debug part dealing with last strip.
C       VX 2.1  JMB     16.6.87 Rectangular images
C       VX 2.2  RH      12.9.87 Extra guide points > 2*IKNOTS per axis.
C       VX 2.3  RH      04.1.88 Extra guide points > 3*IKNOTS per axis.
C       VX 2.4  RH      29.2.88 NDATA=50000
C       VX 2.5  RH      07.4.88 ISIZEX=7500
C       VX 3.0  RH      20.8.88 Efficient guide points &  non-spline option.
C                               renamed CCUNBENDC, old programs retained.
C       VX 3.1  RH     13.11.88 Optimised ISIZEX=6000,ISTEPMAX=60,IDEEP=160
C       VX 3.2  RH     18.12.89 Minor changes, mainly cosmetic.
C       VX 3.3  RH     16.5.90  IMAXCOR limit on X removed.
C       VX 3.4  RH     30.5.90  Fudge on XAPPLY check.
C       VX 3.5  RH     15.6.90  REAL*8 for SUMOUT
C       VX 4.0  RH     03.1.92  Convert for UNIX on Alliant    
C       VX 4.1  RH     23.8.93  insert check for DMIN,DMAX,DMEAN sensible.
C       VX 4.2  RH     07.9.94  increased NDIMX,NDIMY for smoother unbending
C       VX 4.3  RH     09.9.94  improved efficiency for FILLEMPTIES
C       VX 4.4  RH     13.9.94  taperedge option for single molecule work
C                               renames CCUNBENDD, extra input parameters
C       VX 4.5  RH     27.1.95  change plot scale-factor, debug fillempties
C       VX 4.6  RH     29.4.95  add date to plots
C       VX 4.7  RH     25.7.95  ENCODE debug for Alpha
C       VX 4.8  RH     9.10.95  increase dimensions to 8000
C       VX 4.9  RH     23.3.96  optional O/P unbend correction table CCUNBENDE
C       VX 5.0  RH     22.4.96  add initialization and reorder IF statement
C       VX 5.1  RH     11.3.98  change distortion plot to 10x error
C       VX 6.0  RH     23.8.00  convert to plot2000 direct postscript output
C               TSH    13.6.01  P2K_FONT needed string terminator
C
C       MODIFIED  16.6.87  TO BE CORRECT FOR RECTANGULAR IMAGES UP TO SIZE
C       7200 X 9600
C       MODIFIED  14.12.85   TO GREATLY INCREASE THE NUMBER OF KNOTS ALLOWED IN
C            ONLY IN THE DIRECTION PERPENDICULAR TO THE TILT AXIS.
C            WORKS IN CONJUNCTION WITH CCORSERCH OR PROFSERCH.
C            THE PROGRAM UNBENDS THE CRYSTAL USING OUTPUT
C            FROM THE CROSS-CORRELATION PEAK SEARCH PROGRAM, CCORSERCH OR
C            ITS LATER VARIANTS PROFSERCH AND QUADSERCH(0, 1, and 2)
C
C  CONTROL DATA :-
C       CARD 1 : FILE NAME OF INPUT IMAGE FILE.(ONLY HEADER READ IF IOUT=0).
C       CARD 2 : ITYPE,IOUT,IMAXCOR,ISTEP,LTAPER,RTAPER,LTABOUT
CHEN    CARD 3 : IKNOTX,IKNOTY,EPS,FACTOR,TLTAXIS
C       CARD 3 : IKNOTX,IKNOTY,EPS,FACTOR,TLTAXIS,RMAG,LCOLOR
C       CARD 4 : PLOT TITLE FOR DISTORTION CORRECTION DISPLAY.
C      and if IOUT = 1,
C       CARD 5 : FULL FILE NAME FOR OUTPUT OF CORRECTED IMAGE.
C       CARD 6 : TITLE TO BE ADDED TO CORRECTED IMAGE TITLE RECORD.
C
C   INPUT FILES:
C         CCORDATA      - FILE OF PARAMETERS CONTAINING DETAILS OF
C                          THE CCORSERCH RUN AND 
C                       - RAW LIST OF CORRELATION PEAK POSITION AND HEIGHTS
C                          PRODUCED BY CCORSERCH.
C         PIX(NAME)     - ORIGINAL DENSITOMETER RAW IMAGE FILE.
C
C   OUTPUT FILES:
C         PIXOUT(NAME1) - THE UNBENT IMAGE FILE, FULLY CORRECTED FOR
C                          THE SMOOTHED DISTORTION CORRECTIONS.
C         TABLEOUT      - The unbending table as used inside the program - may
C                          be useful for creation of fixed distortion table,
C                          for example, to correct fibre optic distortion in
C                          in another program, e.g. pickprofa.for
C   OPTIONS ARE:
C                         INTERPOLATION OF IMAGE DISTORTION TO PRODUCE
C         ITYPE = 0 ----- USES CORRELATION PEAKS EXACTLY WITHOUT SMOOTHING.
C               = 1 ----- SAME USING BICUBIC SPLINE FITTING (NAGLIB).
C
C         IOUT  = 0 ----- NO OUTPUT, DIAGNOSTICS ONLY.
C               = 1 ----- CORRECTED IMAGE OUTPUT USING THE SMOOTHED PARAMETERS
C                          AS ABOVE.
C         ISTEP --------- SIZE OF GRID ON WHICH THE INTERPOLATED VECTORS ARE
C                          CALCULATED.
C         IMAXCOR ------- SIZE OF THE MAXIMUM ALLOWED CORRECTION, TO ENABLE THE
C                          SIZE OF EACH STRIP READ INTO CORE TO BE CALCULATED -
C                          (ISTEP + 2*IMAXCOR).
C                          CORRECTIONS LARGER THAN IMAXCOR ARE REDUCED.
C                          A FURTHER LIMIT THAT IMAXCOR IS NOT GREATER THAN 
C                          ISTEP WAS INTRODUCED IN Nov-88 TO KEEP THE PROGRAM
C                          SIMPLE. BIGGER VALUES WILL NEED A MAJOR REWRITE.
C         LTAPER -------- T or F (logical*1) for application of a taperedge
C         RTAPER -------- Radius to be used from centre of ISTEP boxes.
C         LTABOUT ------- T or F (logical*1) for output of unbending table
C         IKNOTX -------- NUMBER OF KNOTS TO BE DISTRIBUTED ACROSS IMAGE IN
C                          BICUBIC SPLINE MODE OF SMOOTHING --  IKNOTX DESCRIBES
C                          DISTORTIONS PARALLEL TO THE TILTAXIS.
C         IKNOTY -------- NUMBER OF KNOTS -- IKNOTY DESCRIBES DISTORTIONS
C                          PERPENDICULAR TO THE TILTAXIS -- THIS WILL NEED TO
C                          BE SEVERAL TIMES HIGHER IF THE IMAGE IS OF A HIGHLY
C                          TILTED SPECIMEN.
C         EPS ----------- THRESHOLD FOR DETERMINATION OF RANK OF BICUBIC
C                          SPLINE FITTING MATRIX. TRY 0.00001 --- OTHERWISE
C                          SEE WRITE-UP FOR NAGLIB E02DAF SUBROUTINE.
C         THRESH -------- THRESHOLD OF CROSS-CORRELATION PEAK HEIGHT,
C                          CALCULATED AS;
C                    DENMAX (READ FROM CCORDATA) * FACTOR (READ FROM UNIT 5),
C                        , BELOW WHICH THE PEAK IS NOT USED.
C         TLTAXIS ------- DIRECTION OF TILTAXIS RELATIVE TO NORMAL X-Y AXES OF
C                          IMAGE.
C         RMAG ---------- Magnification for line length in distortion plots
C         CCOLOR ---------y or n (character*1) True for color output in PS plotfiles
C
C*******************************************************************************
C
C  DIMENSION STATEMENTS INCLUDE THE FOLLOWING PARAMETERS :-
C       NDATA  - MAXIMUM NUMBER OF CORRELATION PEAKS ABOVE THRESH.
C       NDIMX   - MAXIMUM NUMBER OF BLOCKS IN WHICH DISTORTION CORRECTION
C                IS CALCULATED. (NUMBER OF BLOCKS = NXYZ(1)/ISTEP),
C                BOTH NXYZ(1) AND ISTEP ARE INPUT PARAMETERS.
C       NDIMY  - NXYZ(2)/ISTEP
C       NMAXKN - MAXIMUM NUMBER OF KNOTS - USABLE NUMBER (IKNOTS)=(NMAXKN-8).
C       NCMAX  -      (NMAXKN+4)**2  SOMETHING TO DO WITH THE KNOT SUBROUTINE.
C       NPOINT - NDATA+(NMAXKN+1)**2       DITTO.
C       NPOINT ALSO - MCALC+(NMAXKN+1)**2; MCALC=NDIMX*NDIMY; MCALC SHOULD NOT 
C               BE LARGER THAN NDATA; VALUE CAN BE ADJUSTED WITH ISTEP
C       ISIZEX  - MAXIMUM IMAGE SIZE IN X-DIMENSION
C       IDEEP  - IMAGE IS READ IN AND CORRECTED IN SLICES OF DEPTH UP TO IDEEP.
C                SMALLER IDEEP RESULTS IN MORE I/O, LARGER IN MORE PAGE FAULTS.
C              - IDEEP MUST BE (2*IMAXCOR + ISTEP); BOTH ON DATA CARD 2
C       NWSPCE - WORKSPACE FOR KNOT SUBROUTINE -- SEE NAGLIB WRITE-UP.
C                  MUST BE .GT. 4+3*(NMAXKN+4)+2*NCMAX*(6+3*(NMAXKN+4))
C                  -- i.e. increases proportional to NMAXKN**3 !!
C       ISTEPMAX - MAXIMUM VALUE OF ISTEP ABLE TO BE USED - CONTROLS DIMENSION
C                  OF PICOUT.
C
      PARAMETER (NDATA=250000)
      PARAMETER (NMAXKN=80)
C-----PARAMETER (NPOINT=NDATA+6561)
      PARAMETER (NPOINT=256561)
      PARAMETER (NCMAX=7056)
      PARAMETER (NWSPCE=410000000)
      PARAMETER (ISIZEX=20100)
      PARAMETER (ISTEPMAX=120)
      PARAMETER (IDEEP=230)
      PARAMETER (NDIMX=2000)
      PARAMETER (NDIMY=2000)
C
C               see subroutine - fillempties, calctaper - parameters also.
      DIMENSION TITLE(20),TITLEIN(20),TITLEERR(20),DATANAME(20)
      REAL*8 SUMOUT
      REAL XPOS(NDATA),YPOS(NDATA),DX(NDATA),DY(NDATA),W(NDATA)
      REAL DXAV(NDIMX,NDIMY),DYAV(NDIMX,NDIMY)        ! Bins for guide point
      REAL NAV(NDIMX,NDIMY),WAV(NDIMX,NDIMY)          ! and linear interpol.
      REAL XOUT(NDIMX*NDIMY),YOUT(NDIMX*NDIMY)
      REAL XCORR(NDIMX*NDIMY),YCORR(NDIMX*NDIMY),WS(NWSPCE)
      REAL LAMBDA(NMAXKN),MU(NMAXKN),DL(NCMAX),CX(NCMAX),CY(NCMAX)
      INTEGER POINT(NPOINT),PX,PY,NCREAL,RANK,IFAIL,J
      REAL MAXCOR, TAPER(ISTEPMAX,ISTEPMAX),TAPER1(ISTEPMAX,ISTEPMAX)
      REAL PICIN(ISIZEX,IDEEP),PICOUT(ISIZEX,ISTEPMAX)
      LOGICAL LTAPER, LTABOUT, LCOLOR
      CHARACTER*1 CCOLOR
      CHARACTER DAT*24
C
      DIMENSION NXYZ(3),MXYZ(3),NXYZ1(3),NXYZST(3)
      CHARACTER*80 NAME,NAME1
      EQUIVALENCE (PICIN(1,1),WS(1))
      EQUIVALENCE (PICOUT(1,1),WS(1+ISIZEX*IDEEP))
      CHARACTER*80 TMPTITLEIN,TMPTITLEERR
      EQUIVALENCE (TMPTITLEIN,TITLEIN),(TMPTITLEERR,TITLEERR)
      DATA NXYZST/3*0/
      DATA NGUIDE/20/         ! MAXIMUM EXTRA GUIDE POINTS = NGUIDE**2.
      DATA TMPTITLEIN/' PLOT OF INPUT DATA ABOVE
     .   THRESH in CCUNBEND K'/
      DATA TMPTITLEERR/' PLOT OF FITTING ERROR AT
     .   INPUT DATA POINTS in CCUNBEND2K'/
C
      XCOORD(I,J)=A1*I+B1*J+IC
      YCOORD(I,J)=A2*I+B2*J+IR
C*** initialization added by jms 06.03.96
      do j=1,ndimy
        do i=1,ndimy
          WAV(i,j) = 0.
          NAV(i,j) = 0
          DXAV(i,j) = 0.
          DYAV(i,j) = 0.
        end do
      end do
      do j=1,istepmax
        do i=1,istepmax
          taper(i,j) = 0.
        end do
      end do
C
C  INPUT OF LIST OF ALL IMAGE PARAMETERS AS IN PRODUCED BY CCORSERCH
      WRITE(6,13211)
13211 FORMAT(/,/,' CCUNBENDK VX 6.0(23.8.00) - program to reinterpolate',
     . ' an image of a crystal on a straight lattice',/,/)
      CALL CCPDPN(3,'CCORDATA','READONLY','F',0,0)
      READ(3,13210)DATANAME
13210 FORMAT (20A4)
      WRITE(6,3210)DATANAME
      READ(3,13210)DATANAME
      WRITE(6,3210)DATANAME
      READ(3,13210)DATANAME
      WRITE(6,3210)DATANAME
C
      READ(3,*) NC,NR,IC,IR,A1,A2,B1,B2,MINA,MAXA,MINB,MAXB
      READ(3,*)DENMAX
      WRITE(6,9001)NC,NR,IC,IR,A1,A2,B1,B2,MINA,MAXA,MINB,MAXB,DENMAX
9001  FORMAT('$NC,NR ',2I5/
     .  '$IC,IR ',2I5/
     .  '$A1,A2,B1,B2 ',4F10.4/
     .  '$LAST XCOR X,Y ',4I5/
     .  '$MAX DENSITY IN CORRELATION MAP',F12.1)
C READ IN NAME OF IMAGE FILE TO BE USED.
      WRITE(6,3211)
      READ(5,163)NAME
      WRITE(6,3212)NAME
C
      CALL FDATE(DAT)
      WRITE(6,'('' Date is '',A20)') DAT(5:24)
C
      call shorten(NAME,k)
      if(k.gt.60)k=60
      WRITE(TMPTITLEIN(1:80),'(80X)')
      WRITE(TMPTITLEIN(1:60),'(A)')NAME(1:k)
      WRITE(TMPTITLEIN(61:80),'(A)')DAT(5:24)
C
C  INPUT OF ALL CONTROL DATA REQUIRED TO RUN THIS PROGRAM.
      WRITE(6,160)
      READ(5,*)ITYPE,IOUT,IMAXCOR,ISTEP,LTAPER,RTAPER,LTABOUT
      IF(IMAXCOR.GT.ISTEP) THEN
        WRITE(6,159) IMAXCOR,ISTEP
        IMAXCOR=ISTEP
      ENDIF
      WRITE(6,161)ITYPE,IOUT,IMAXCOR,ISTEP,LTAPER,RTAPER,LTABOUT
      WRITE(6,169)
      READ(5,*) IKNOTX,IKNOTY,EPS,FACTOR,TLTAXIS,RMAG,CCOLOR
      if(CCOLOR.eq."y")then
        LCOLOR=.true.
      else
        LCOLOR=.false.
      endif
      MAXCOR=IMAXCOR-0.0001
      THRESH=FACTOR*DENMAX
      WRITE(6,168) IKNOTX,IKNOTY,EPS,FACTOR,THRESH,TLTAXIS
C
C  READ IMAGE HEADER TO CHECK THAT REQUESTED PARAMETERS AND IMAGE SIZE ARE NOT
C  TOO LARGE FOR DECLARED MAXIMUM PROGRAM DIMENSIONS (IN PARAMETER STATEMENT).
      CALL IMOPEN(1,NAME,'RO')
      CALL IRDHDR(1,NXYZ,MXYZ,MODE,DMIN,DMAX,DMEAN)
CHENN>
C      IF(DMEAN.LT.0.0.OR.DMEAN.GT.5000.0
C     .   .OR.DMIN.LT.1.0.OR.DMAX.LT.1.0) THEN
C      IF(DMEAN.LT.0.0.OR.DMEAN.GT.50000.0
C     . .OR.DMIN.LT.-10000.0.OR.DMAX.LT.1.0) THEN
CHENN<
C        WRITE(6,7500)
C7500    FORMAT(' ',/,':WARNING: Image values for DMIN,DMAX or',
C     .    ' DMEAN are:')
C        write(6,'('':DMIN = '',G15.3,'', DMAX = '',G15.3,
C     .            '', DMEAN= '',G15.3)')DMIN,DMAX,DMEAN
CCHEN>
C       STOP
CHEN<
C      ENDIF
C
C  CHECK NOW THAT PROGRAM DIMENSIONS ARE ADEQUATE.
      IF(ISIZEX.LT.NXYZ(1)) THEN
        IDIMOUT=1
        GO TO 199
      END IF
      IF(ISTEP.GT.ISTEPMAX) THEN
        IDIMOUT=2
        GO TO 199
      END IF
      IF(NDIMX.LT.(NXYZ(1)/ISTEP)) THEN
        IDIMOUT=3
        GO TO 199
      END IF
      IF(NDIMY.LT.(NXYZ(2)/ISTEP)) THEN
        IDIMOUT=4
        GO TO 199
      END IF
C
C  CALCULATE COMPRESSION OF CORRELATION MAP WITH RESPECT TO REAL IMAGE
C   (eg 1,2 OR 3)
C
      ANC=NC
      ANR=NR
      COMPRSSX=NXYZ(1)/ANC
      COMPRSSY=NXYZ(2)/ANR
      WRITE(6,15000)COMPRSSX,COMPRSSY
15000 FORMAT(/,' RATIO ( NUMBER OF SAMPLES IN CORRELATION MAP)',/,
     . '          -----------------------------------',/,
     . '          ( NUMBER OF SAMPLES IN REAL IMAGE )',/,/,
     . '     COMPRESS IN X DIRECTION = ',F8.5,/,
     . '                 Y DIRECTION = ',F8.5,/,/)
      ICOMPRSSX=NINT(COMPRSSX)
      ICOMPRSSY=NINT(COMPRSSY)
      WRITE(6,14999)ICOMPRSSX,ICOMPRSSY
14999 FORMAT(' ICOMPRSSX=',I5,' ICOMPRSSY=',I5,/)
      IF(ITYPE.EQ.0) GO TO 250
      IF(IKNOTX.GT.(NMAXKN-8)) GO TO 197
      IF(IKNOTY.GT.(NMAXKN-8)) GO TO 197
      IF(NPOINT.LT.(NDATA+(NMAXKN+1)**2)) THEN
        IDIMOUT=5
        GO TO 199
      END IF
      IF(NCMAX.LT.(NMAXKN+4)**2) THEN
        IDIMOUT=6
        GO TO 199
      END IF
      NTEST=4+3*(NMAXKN+4)+2*NCMAX*(6+3*(NMAXKN+4))
      IF(NWSPCE.LT.NTEST) THEN
        IDIMOUT=7
        GO TO 199
      END IF
250   CONTINUE
C
C  READ IN CCORDATA DATA AND STORE NON-ZERO CORRELATION POSITIONS
      MDATA=0
      NLOST=0
      NOUTSIDE=0
C
      XMIN = NXYZ(1)
      XMAX = 0.0
      YMIN = NXYZ(2)
      YMAX = 0.0
      DISTMAX = 0.0
C      DO 200 I=MINA,MAXA
C        DO 200 J=MINB,MAXB
      DO 200 Itmp=MINA,MAXA
        DO 200 Jtmp=MINB,MAXB
C          READ(3,*) XCOOR,YCOOR,PEAK
          READ(3,*,END=203) I,J,XCOOR,YCOOR,PEAK
          IF(XCOOR.EQ.0.0) GO TO 200
C          write(*,'('' X/YCOOR = '',2F10.1,''  PEAK = '',F10.3)')
C     .      XCOOR,YCOOR,PEAK
          IF(PEAK.LT.THRESH) THEN
            NLOST=NLOST+1
            GO TO 200
          ENDIF
          MDATA=MDATA+1
          IF(MDATA.GT.NDATA) THEN
            WRITE(6,104)MDATA
            STOP
          ENDIF
          XPOS(MDATA)=XCOORD(I,J)*ICOMPRSSX
          YPOS(MDATA)=YCOORD(I,J)*ICOMPRSSY
          IF(XPOS(MDATA).LT.1.0.OR.XPOS(MDATA).GT.NXYZ(1).OR.
     .      YPOS(MDATA).LT.1.0.OR.YPOS(MDATA).GT.NXYZ(2)) THEN
CHENN>
C           WRITE(6,276) XPOS(MDATA),YPOS(MDATA),I,J
C276        FORMAT(' !!!!!!!  -- point outside image ignored --',
C     .       'XPOS(I),YPOS(I),I,J,MDATA= ',2F10.2,3I6)
            WRITE(6,276) XPOS(MDATA),YPOS(MDATA),I,J,MDATA,XCOOR,YCOOR
276         FORMAT(' !!!!!!!  -- point outside image ignored --',
     .        'XPOS(I),YPOS(I),I,J,MDATA= ',2F10.2,3I6,2F10.2)
CHENN<
            NOUTSIDE=NOUTSIDE+1
            MDATA=MDATA-1
            GO TO 200
          ENDIF
C ----------------------------------------------------------------
          XMIN=AMIN1(XPOS(MDATA),XMIN)
          XMAX=AMAX1(XPOS(MDATA),XMAX)
          YMIN=AMIN1(YPOS(MDATA),YMIN)
          YMAX=AMAX1(YPOS(MDATA),YMAX)
          DX(MDATA)=(XCOOR-XCOORD(I,J))*ICOMPRSSX
          DY(MDATA)=(YCOOR-YCOORD(I,J))*ICOMPRSSY
          W(MDATA)=PEAK         ! WEIGHT PROPNL TO CORREL PEAK HEIGHT.
          IBINX=1+XPOS(MDATA)/ISTEP
          IBINY=1+YPOS(MDATA)/ISTEP
          WAV(IBINX,IBINY) = WAV(IBINX,IBINY)+W(MDATA)
          DXAV(IBINX,IBINY)= DXAV(IBINX,IBINY)+DX(MDATA)
          DYAV(IBINX,IBINY)= DYAV(IBINX,IBINY)+DY(MDATA)
          NAV(IBINX,IBINY) = NAV(IBINX,IBINY) + 1
          DISTCORR = SQRT(DX(MDATA)**2 + DY(MDATA)**2)
          IF(DISTCORR.GT.DISTMAX) DISTMAX = DISTCORR
200   CONTINUE
203   continue
C
C  CALCULATE AT HOW MANY POSITIONS SMOOTHED DISTORTION WILL BE DETERMINED.
      MXDIM=(NXYZ(1)+ISTEP-1)/ISTEP
      MYDIM=(NXYZ(2)+ISTEP-1)/ISTEP
      MCALC=MXDIM*MYDIM
C
C  NOW FIND THE MEAN CORRECTION AND WEIGHT IN EACH OF THE BIN BOXES
      DO 260 I=1,MXDIM
        DO 260 J=1,MYDIM
          IF(NAV(I,J).EQ.0) GO TO 260
          WAV(I,J) = WAV(I,J)/NAV(I,J)
          DXAV(I,J)= DXAV(I,J)/NAV(I,J)
          DYAV(I,J)= DYAV(I,J)/NAV(I,J)
260   CONTINUE
C
      WRITE(6,261)
261   FORMAT('  STARTING SUBROUTINE FILLEMPTIES')
      CALL FILLEMPTIES(MXDIM,MYDIM,DXAV,DYAV,WAV,NAV)
C
      WRITE(6,100)MDATA,NLOST,DISTMAX,NOUTSIDE
      MDATAOBS=MDATA
      CALL PLOTCORR(TITLEIN,MDATA,XPOS,YPOS,DX,DY,NXYZ,100,RMAG,LCOLOR)
      IF(ITYPE.EQ.1)CALL ROTATEXY(MDATA,DX,DY,TLTAXIS)  ! BICUBIC ONLY
C
C  ABOVE PART OF PROGRAM HAS CALCULATED THE VALUES OF THE DISTORTION
C  CORRECTIONS AT EACH OF A NUMBER OF POINTS IN THE IMAGE. BELOW A FEW 
C  POINTS ARE ADDED TO THE IMAGE SO THAT A TRUE BOUNDARY IS CREATED FOR
C  INTERPOLATION AND SO THAT THERE ARE NO EMPTY SPACES IN THE IMAGE.
C  RESTRICT NUMBER OF EXTRA POINTS TO NGUIDE**2 FOR SPEED.
C
        MEDGE=SQRT(FLOAT(MDATA))
        IF(ITYPE.EQ.1) THEN
                MEDGE=MAX(MEDGE,NGUIDE,3*IKNOTX,3*IKNOTY)
        ELSE 
                MEDGE=MAX(MEDGE,NGUIDE)
        ENDIF
          NGUIDMIN=MAX(1,(MEDGE**2*ISTEP**2/(NXYZ(1)*NXYZ(2))))
          WRITE(6,108) NGUIDMIN ! unlikely to be different from 1,in most cases.
        XINT=FLOAT(NXYZ(1)-1)/(MEDGE-1)
        YINT=FLOAT(NXYZ(2)-1)/(MEDGE-1)
        DO 270 I=1,MEDGE
                XTEST=1+(I-1)*XINT
                IF(I.EQ.MEDGE) XTEST=NXYZ(1)    ! ensures true edge is reached
        DO 270 J=1,MEDGE
                YTEST=1+(J-1)*YINT
                IF(J.EQ.MEDGE) YTEST=NXYZ(2)    ! ensures true edge is reached
                IBINX=1+XTEST/ISTEP
                IBINY=1+YTEST/ISTEP
C
C          GUIDE POINTS ROUND EDGE WILL ALWAYS BE ADDED TO
C          ENSURE THAT CUBIC SPLINE ROUTINES ARE OK.
                IF(I.EQ.1.OR.I.EQ.MEDGE.OR.J.EQ.1.OR.J.EQ.MEDGE)GOTO 275
                IF(NAV(IBINX,IBINY).GE.NGUIDMIN) GO TO 270  ! i.e. no guide pt.
275     CONTINUE
        MDATA=MDATA+1
      IF(MDATA.GT.NDATA) THEN
        WRITE(6,104)MDATA
        STOP
      ENDIF
        XPOS(MDATA)=XTEST       ! Adding single guide point per ISTEPxISTEP box
        YPOS(MDATA)=YTEST       !   "     "
C  MAKE DISTORTION AT EDGE AND FILL-IN POINTS SAME AS AVERAGE OF BINNED INPUT
C  DATA POINTS.
        DX(MDATA) = DXAV(IBINX,IBINY)
        DY(MDATA) = DYAV(IBINX,IBINY)
        W(MDATA)  = WAV(IBINX,IBINY)/4.0        ! but lower weights
270     CONTINUE
      MADDED=MDATA-MDATAOBS
      WRITE(6,105)MDATA,MADDED
C
C  CALCULATE POSITIONS AT WHICH SMOOTHED DISTORTION WILL BE DETERMINED.
        XMAX=NXYZ(1)
        YMAX=NXYZ(2)
        ICALC=0
      DO 400 J=1,MYDIM
      DO 400 I=1,MXDIM
        ICALC=ICALC+1
        XOUT(ICALC)=ISTEP*(I-0.5)
        IF(XOUT(ICALC).GT.XMAX) XOUT(ICALC)=XMAX
        YOUT(ICALC)=ISTEP*(J-0.5)
        IF(YOUT(ICALC).GT.YMAX) YOUT(ICALC)=YMAX
400   CONTINUE
      WRITE(6,107)MCALC,ICALC
C
C  this section for debugging a problem WITH POSITIONS OF KNOTS.
C
        XDATMIN=1000
        XDATMAX=0
        YDATMIN=1000
        YDATMAX=0
        DO 405 I=1,MDATA
        XDATMIN=AMIN1(XDATMIN,XPOS(I))
        XDATMAX=AMAX1(XDATMAX,XPOS(I))
        YDATMIN=AMIN1(YDATMIN,YPOS(I))
405     YDATMAX=AMAX1(YDATMAX,YPOS(I))
        XOUTMIN=1000
        XOUTMAX=0
        YOUTMIN=1000
        YOUTMAX=0
        DO 406 I=1,MCALC
        XOUTMIN=AMIN1(XOUTMIN,XOUT(I))
        XOUTMAX=AMAX1(XOUTMAX,XOUT(I))
        YOUTMIN=AMIN1(YOUTMIN,YOUT(I))
406     YOUTMAX=AMAX1(YOUTMAX,YOUT(I))
      WRITE(6,407) XDATMIN,XDATMAX,YDATMIN,YDATMAX,
     . XOUTMIN,XOUTMAX,YOUTMIN,YOUTMAX
407     FORMAT( ' MIN & MAX X & Y on input+guide data points =',4F10.3/
     . ' min & max X & Y at centres of ISTEP boxes  =',4F10.3)
      IF(XOUTMAX.GT.XDATMAX.OR.YOUTMAX.GT.YDATMAX) THEN
                WRITE(6,4060)XOUTMAX,XDATMAX,YOUTMAX,YDATMAX
4060            FORMAT(///' ERROR IN MAXIMUM CALC GUIDE PT',4F15.5)
                STOP
      ENDIF
C   END OF DEBUG BIT
C
C
C
C
C  SIMPLE CORRECTION BY THE AMOUNT OF THE LOCAL CORRELATION PEAK.
      IF(ITYPE.EQ.0) THEN
        DO 410 I=1,MCALC
          IBINX=1+XOUT(I)/ISTEP
          IBINY=1+YOUT(I)/ISTEP
          XCORR(I)=DXAV(IBINX,IBINY)
          YCORR(I)=DYAV(IBINX,IBINY)
410     CONTINUE
      ENDIF
C
C  MORE COMPLICATED CORRECTION WITH BICUBIC SPLINE SMOOTHING.
      IF(ITYPE.EQ.1) THEN                       ! ENDS ABOUT 99 LINES DOWN.
        IF(NPOINT.LT.(MCALC+(NMAXKN+1)**2)) THEN
          IDIMOUT=8
          GO TO 199
        END IF
C
C  CALCULATE POSITIONS OF BICUBIC SPLINE KNOTS FIRST FROM IKNOTX AND
C   THEN LATER (SEE BELOW) FOR IKNOTY.
C
        PX=8+IKNOTX
        PY=8+IKNOTX
        DO 300 I=1,IKNOTX
          LAMBDA(I+4)=(NXYZ(1)/(IKNOTX+1))*I
300     MU(I+4)=(NXYZ(2)/(IKNOTX+1))*I
        NCREAL=(PX-4)*(PY-4)
        IFAIL=0
        NADRES=(PX-7)*(PY-7)
C  NOW SORT SPOTS
        CALL E02ZAF(PX,PY,LAMBDA,MU,MDATA,XPOS,YPOS,POINT,NPOINT,WS,
     .  NADRES,IFAIL)
        IF(IFAIL.NE.0) THEN
          WRITE(6,101)IFAIL
          STOP
        ENDIF
C  NOW DO SPLINE FIT FIRST ON X-ERROR (IE. PARALLEL TO TILTAXIS) == DX
        CALL E02DAF(MDATA,PX,PY,XPOS,YPOS,DX,W,LAMBDA,MU,POINT,NPOINT,
     .   DL,CX,NCREAL,WS,NWSPCE,EPS,SIGMA,RANK,IFAIL)
        IF(IFAIL.NE.0) THEN
          WRITE(6,102)IFAIL
          STOP
        ENDIF
        WRITE(6,150)SIGMA 
        WRITE(6,151)RANK,NCREAL,EPS
C
C  NOW CALCULATE SMOOTHED SPLINE FIT AT FINE ENOUGH INTERVALS FOR OUTPUT
C  PLOTTING AND FOR IMAGE CORRECTION - AS CALCULATED ABOVE.
C
        CALL E02ZAF(PX,PY,LAMBDA,MU,MCALC,XOUT,YOUT,POINT,NPOINT,WS,
     .    NADRES,IFAIL)
        IF(IFAIL.NE.0) THEN
          WRITE(6,1101)IFAIL
          STOP
        ENDIF
C  CREATION OF XCORR
        WRITE(6,408)LAMBDA(4),LAMBDA(PX-3),MU(4),MU(PY-3)
408     FORMAT(' LAMBDA(4,PX-3) & MU(",")',4F8.1)
        WRITE(6,409) (J,LAMBDA(J),MU(J),J=1,PX)
409     FORMAT(' ALL LAMBDAS & MUS, N   LAMBDA   MU'/90(15X,I5,2F8.1/))
        CALL E02DBF(MCALC,PX,PY,XOUT,YOUT,XCORR,LAMBDA,MU,POINT,NPOINT,
     .    CX,NCREAL,IFAIL)
        IF(IFAIL.NE.0) THEN
          WRITE(6,1102)IFAIL
          STOP
        ENDIF
C
C  NOW CALCULATE POSITIONS OF BICUBIC SPLINE KNOTS FOR IKNOTY.
        PX=8+IKNOTY
        PY=8+IKNOTY
        DO 301 I=1,IKNOTY
          LAMBDA(I+4)=(NXYZ(1)/(IKNOTY+1))*I
301     MU(I+4)=(NXYZ(2)/(IKNOTY+1))*I
        NCREAL=(PX-4)*(PY-4)
        IFAIL=0
        NADRES=(PX-7)*(PY-7)
C  NOW SORT SPOTS AS FOR X ABOVE
        CALL E02ZAF(PX,PY,LAMBDA,MU,MDATA,XPOS,YPOS,POINT,NPOINT,WS,
     .    NADRES,IFAIL)
        IF(IFAIL.NE.0) THEN
          WRITE(6,101)IFAIL
          STOP
        ENDIF
C  NOW ON Y-ERROR DY
        CALL E02DAF(MDATA,PX,PY,XPOS,YPOS,DY,W,LAMBDA,MU,POINT,NPOINT,
     .     DL,CY,NCREAL,WS,NWSPCE,EPS,SIGMA,RANK,IFAIL)
        IF(IFAIL.NE.0) THEN
          WRITE(6,103)IFAIL
          STOP
        ENDIF
        WRITE(6,150)SIGMA
        WRITE(6,151)RANK,NCREAL,EPS
C
C  NOW CALCULATE SMOOTHED SPLINE FIT AT FINE ENOUGH INTERVALS FOR OUTPUT
C  PLOTTING AND FOR IMAGE CORRECTION - AGAIN AS ABOVE.
C
        CALL E02ZAF(PX,PY,LAMBDA,MU,MCALC,XOUT,YOUT,POINT,NPOINT,WS,
     .    NADRES,IFAIL)
        IF(IFAIL.NE.0) THEN
          WRITE(6,1101)IFAIL
          STOP
        ENDIF
C  NOW CREATION OF YCORR AS FOR XCORR ABOVE
        CALL E02DBF(MCALC,PX,PY,XOUT,YOUT,YCORR,LAMBDA,MU,POINT,NPOINT,
     .    CY,NCREAL,IFAIL)
        IF(IFAIL.NE.0) THEN
          WRITE(6,1103)IFAIL
          STOP
        ENDIF
C
      ENDIF   !    END OF CUBIC SPLINE FITTING (ITYPE=1).
C
C      BELOW IS PLOTTING OF THE FITTED DEVIATIONS AT ALL CALCULATED POSNS.
C
1107  FORMAT(' TITLE FOR PLOT OF SMOOTHED CORRECTIONS ?'/'$      :-')
1108  FORMAT(20A4)
1109  FORMAT(' TITLE WAS  ',20A4)
      WRITE(6,1107)
      READ(5,1108)  TITLE
      WRITE(6,1109) TITLE
      IF(ITYPE.EQ.1) THEN               ! BICUBIC ONLY
        CALL ROTATEXY(MDATAOBS,DX,DY,-TLTAXIS)
        CALL ROTATEXY(MCALC,XCORR,YCORR,-TLTAXIS)
      ENDIF
      CALL PLOTCORR(TITLE,MCALC,XOUT,YOUT,XCORR,YCORR,NXYZ,1,RMAG,LCOLOR)
      IF(LTABOUT) 
     . CALL WRITETABLE(TITLE,ISTEP,MXDIM,MYDIM,NXYZ,XCORR,YCORR)

C
C  NOW OUTPUT OF FITTING ERROR AT INPUT DATA POINTS
C
C  CALCULATE NEAREST POINT IN SMOOTHED FUNCTION
C USE DX, DY ARRAYS TO CALCULATE THE RESIDUAL FITTING ERROR FOR PLOTTING.
      STEPROOT2=ISTEP/SQRT(1.999)
      WRITE(6,1106)STEPROOT2
      NOTCORRFUL=0
      DO 1104 I=1,MDATAOBS
        K= 1 + (XPOS(I)/ISTEP) + MXDIM*INT(YPOS(I)/ISTEP)
        DIST = (XOUT(K)-XPOS(I))**2 + (YOUT(K)-YPOS(I))**2
        DIST = SQRT(DIST)
        IF(DIST.GT.STEPROOT2) THEN              ! CHECK CALC IS CORRECT
          WRITE(6,1105)I,XPOS(I),YPOS(I),K,XOUT(K),YOUT(K)
        ENDIF
        XAPPLY=XCORR(K)
        YAPPLY=YCORR(K)
        IF(ABS(YAPPLY).GT.MAXCOR) YAPPLY=SIGN(MAXCOR,YAPPLY)
        DX(I)=DX(I)-XAPPLY
        DY(I)=DY(I)-YAPPLY
        IF(ABS(YCORR(K)).GT.MAXCOR)
     .   NOTCORRFUL=NOTCORRFUL+1                ! count up truncations
1104  CONTINUE
      if(MDATAOBS.gt.0)then
        PROPNLOST=FLOAT(NOTCORRFUL)/MDATAOBS
      else
        PROPNLOST=100.0
      endif
      if(PROPNLOST.gt.0.1)then
        WRITE(6,1111)IMAXCOR,NOTCORRFUL,PROPNLOST
      else
        WRITE(6,1112)IMAXCOR,NOTCORRFUL,PROPNLOST
      endif
1111  FORMAT(':: THERE WAS A TOTAL OF CORRECTION VECTORS OF LENGTH',/,
     .  ':: GREATER THAN IMAXCOR(=',I4,')  OF',I8,10X,
     .  ' !!!! BEWARE DATA LOST !!!!',/,':: ***** PROPORTION LOST =',F7.3)
1112  FORMAT(' THERE WAS A TOTAL OF CORRECTION VECTORS OF LENGTH',/,
     .  ' GREATER THAN IMAXCOR(=',I4,')  OF',I8,10X,
     .  ' (This is lost data)',/,' ***** PROPORTION LOST =',F7.3)
      CALL PLOTCORR(TITLEERR,MDATAOBS,XPOS,YPOS,DX,DY,NXYZ,-1,RMAG,LCOLOR)
C
      IF(IOUT.EQ.0) STOP
      IF((ISTEP+2*IMAXCOR).GT.IDEEP) GO TO 198
C
CHENN>
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  NOW APPLY THIS SMOOTHED CORRECTION TO THE IMAGE DENSITY VALUES TO PRODUCE
C  A CORRECTED IMAGE.
C
C  TRANSFER OLD HEADER TO OUTPUT FILE.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CHENN<
C
      NXYZ1(1) = NXYZ(1)
      NXYZ1(2) = NXYZ(2)
      NXYZ1(3) = 1
      WRITE(6,162)
      READ(5,163) NAME1
      WRITE(6,164) NAME1
      WRITE(6,165)
      READ(5,166) TITLE
      WRITE(6,167) TITLE
      CALL IMOPEN(4,NAME1,'NEW')
      CALL ITRHDR(4,1)
      CALL IALSIZ(4,NXYZ1,NXYZST)
      CALL IWRHDR(4,TITLE,1,DMIN,DMAX,DMEAN)
C
C  READ IN A STRIP OF SEVERAL LINES   (ISTEP + 2*IMAXCOR)   DEEP.
C
      NJUMPX=0
      NJUMPY=0
      NCHANGE=0
      SUMOUT=0.0
      DMINOUT=DMEAN
      DMAXOUT=DMEAN
      ICALC=0
      NX1=0
      NX2=NXYZ(1)-1
      DO 500 J=1,MYDIM
C  INPUT OF STRIP OF IMAGE...................
C  READ IN REGION (ISTEP + 2*IMAXCOR) DEEP, EXCEPT FOR THE FIRST AND
C   LAST STRIPS WHICH BEGIN AND END AT THE PICTURE EDGE.
        CALL IMPOSN(1,0,0)
        NY1 =  (J-1)*ISTEP -IMAXCOR     ! EXTRA IMAXCOR READ IN PER STRIP
        IF(J.EQ.1) NY1 = 0      ! EXCEPT FOR FIRST STRIP
        NY2 = J*ISTEP + IMAXCOR - 1                     ! SAME AT END OF STRIP
        IF(NY2.GT.NXYZ(2)-1) THEN           ! EXCEPT FOR LAST STRIP(S)
          NY2KEEP=NY2
          NY2 = NXYZ(2)-1   
        END IF
        CALL IRDPAS(1,PICIN,ISIZEX,IDEEP,NX1,NX2,NY1,NY2,*9500)
C
        ITEST=ISTEP+2*IMAXCOR   ! TEST FOR Y OUTSIDE RANGE OF PICIN.
        IF(J.EQ.1) ITEST=ISTEP+IMAXCOR                             !FIRST STRIP
        IF(NY2KEEP.GT.NXYZ(2)-1) ITEST=IMAXCOR+NXYZ(2)-(J-1)*ISTEP !LAST STRIPS
C
C SET PICOUT TO PICIN SO POINTS, WHICH (AFTER DISTORTION CORRECTION)
C ARE OFF THE EDGE OF THE ORIGINAL IMAGE, ARE SET TO ORIGINAL VALUE.
C
        DO 490 IX=1,NXYZ(1)
          DO 490 IY=1,ISTEP
            PICOUT(IX,IY)=DMEAN
490     continue
        IF(LTAPER) CALL CALCTAPER(ISTEP,TAPER,RTAPER)
        DO 493 IS=1,ISTEP
          DO 493 JS=1,ISTEP
            TAPER1(IS,JS)=DMEAN*(1.0-TAPER(IS,JS))
493     continue
C
C  DO BILINEAR INTERPOLATION OF IMAGE DENSITY  :  PICIN --> PICOUT.
C
C
        ierrorcount1=0
        ierrorcount2=0
        DO 480 I=1,MXDIM
          ICALC=ICALC+1
          DO 450 IX=1,ISTEP
            IXTRUE=IX + (I-1)*ISTEP
            IF (IXTRUE.GT.NXYZ(1))GO TO 450
            DO 452 IY=1,ISTEP
              IYTRUE=IY + (J-1)*ISTEP
              IF (IYTRUE.GT.NXYZ(2))GO TO 452
              IYINPIC = IY + IMAXCOR    ! PICIN STARTS AT (IYTRUE-IMAXCOR).
              IF(J.EQ.1) IYINPIC=IY     ! EXCEPT FOR FIRST STRIP
              XAPPLY=XCORR(ICALC)
              YAPPLY=YCORR(ICALC)
              IF(ABS(YAPPLY).GT.MAXCOR) YAPPLY=SIGN(MAXCOR,YAPPLY)
              IXA=XAPPLY+IXTRUE
              IYA=YAPPLY+IYINPIC
              IF((IXA.LT.1).OR.(IXA.GT.(NXYZ(1)-1))) THEN  !TEST FOR X OUTSIDE IMAGE
                IF(I.EQ.1.OR.I.EQ.MXDIM.OR.I.EQ.MXDIM-1)
     .            GO TO 1450 ! PROB OK BUT NO CHECK POSS.
                IF(I.LE.3.OR.I.GE.MXDIM-3)GO TO 1450 !Fudge for removing Xmaxcor
                if(ierrorcount1.lt.30)then
                  WRITE(6,191)IXA,(NXYZ(1)-1),IX,XAPPLY         ! ERROR
                  WRITE(6,193)I,J,IX,IY,IXTRUE,IYTRUE           ! ERROR
                  ierrorcount1=ierrorcount1+1
                else
                  if(ierrorcount1.eq.30)then
                    write(6,'(/,/,/,''ERROR NOT FURTHER REPORTED'',/,/)')
                  endif
                    ierrorcount1=ierrorcount1+1
                endif
                GO TO 1450
              ENDIF
              IF((IYA.LT.1).OR.(IYA.GT.(ITEST-1))) THEN ! ITEST DEFINED ABOVE
                IF(J.EQ.1.OR.NY2KEEP.GT.NXYZ(2)-1) GO TO 1451
C                                               ! OK BUT NO CORRECTION POSSIBLE
                if(ierrorcount2.lt.30)then
                  WRITE(6,192)IYA,ITEST,IDEEP,IY,YAPPLY ! ERROR
                  WRITE(6,193)I,J,IX,IY,IXTRUE,IYTRUE   ! ERROR
                  ierrorcount2=ierrorcount2+1
                else
                  if(ierrorcount2.eq.30)then
                    write(6,'(/,/,/,''ERROR NOT FURTHER REPORTED'',/,/)')
                  endif
                  ierrorcount2=ierrorcount2+1
                endif
                GO TO 1451
              ENDIF
              XDELTA=XAPPLY+IXTRUE-IXA
              YDELTA=YAPPLY+IYINPIC-IYA
              PICOUT(IXTRUE,IY)=(PICIN(IXA,IYA) * (1.-XDELTA)
     .        +  PICIN(IXA+1,IYA) * XDELTA) * (1.0-YDELTA)
     .        +  (PICIN(IXA,IYA+1) * (1.-XDELTA)
     .        +  PICIN(IXA+1,IYA+1) * XDELTA) * YDELTA
              IF(LTAPER) PICOUT(IXTRUE,IY) = 
     .          PICOUT(IXTRUE,IY)*TAPER(IX,IY) + 
     .          TAPER1(IX,IY)
              IF(PICOUT(IXTRUE,IY).LT.DMINOUT)THEN
                IXMINOUT=IXTRUE
                IYMINOUT=IYTRUE
                DMINOUT=PICOUT(IXTRUE,IY)
              ENDIF
              IF(PICOUT(IXTRUE,IY).GT.DMAXOUT)DMAXOUT=PICOUT(IXTRUE,IY)
              NCHANGE=NCHANGE+1
              GO TO 451
1450          NJUMPX=NJUMPX+1
              GO TO 451
1451          NJUMPY=NJUMPY+1
451           SUMOUT=SUMOUT+PICOUT(IXTRUE,IY)
452         CONTINUE
450       CONTINUE
480     CONTINUE
C
C  OUTPUT OF STRIP OF CORRECTED UNBENT IMAGE.
        IF(J.EQ.MYDIM) THEN
          INOUT=NXYZ(2)-(J-1)*ISTEP
        ELSE
          INOUT=ISTEP
        ENDIF
CHEN    WRITE(6,195) J,INOUT
        DO 495 M=1,INOUT
          CALL IWRLIN(4,PICOUT(1,M))
495     CONTINUE
500   CONTINUE
      CALL IMCLOSE(1)
      DMEANOUT=SUMOUT/(NCHANGE+NJUMPX+NJUMPY)
      WRITE(6,20001)DMINOUT,DMAXOUT,DMEANOUT
20001 FORMAT(/' OUTPUT DENSITIES; MIN=',E10.4,' MAX=',E10.4,
     . ' MEAN=',E10.4/)
      WRITE(6,20003)IXMINOUT,IYMINOUT
20003 FORMAT(' COORDS OF MIN DENSITY',2I10/)
      WRITE(6,20002)NCHANGE,NJUMPX,NJUMPY
20002 FORMAT(1X,I10,' ALTERED DENSITIES AND',2I7,' UNCHANGED'/)
      CALL IWRHDR(4,TITLE,-1,DMINOUT,DMAXOUT,DMEANOUT)
      CALL IMCLOSE(4)
C
      STOP
9500    WRITE(6,9012)NX1,NX2,NY1,NY2,ISIZEX,IDEEP
9012    FORMAT(':: Error reading picture file from IRDPAS'/
     . '::   Trying to read a strip from',I6,' to',I6,' in X'/
     . '::                      and from',I6,' to',I6,' in Y'/
     . ':: into array PICIN of dimension',I6,' by',I6)
        STOP
197     WRITE(6,172) NMAXKN-8
        STOP
198     WRITE(6,171)
        STOP
199     WRITE(6,170)IDIMOUT
        STOP
100     FORMAT(' NUMBER OF CROSS-CORRELATION PEAKS READ IN; MDATA',I10/
     . ' NUMBER OF PEAKS BELOW THRESHOLD NOT READ IN; NLOST',I10/
     . ' MAXIMUM DISTORTION CORRECTION READ IN; DISTMAX',F20.1/
     . ' NUMBER OF POINTS NOT READ IN BECAUSE THEY WOULD BE OUTSIDE',
     . ' IMAGE AFTER UNBENDING; NOUTSIDE',I10)
101     FORMAT(':: FAILED KNOT SORTING E02ZAF, IFAIL=',I5)
102     FORMAT(':: CUBIC SPLINE FIT OF X-ERROR FAILED, IFAIL=',I5)
103     FORMAT(':: CUBIC SPLINE FIT OF Y-ERROR FAILED, IFAIL=',I5)
104     FORMAT(':: NDATA PROGRAM DIMENSION TOO SMALL',I10)
105     FORMAT(' TOTAL OF MDATA AFTER ADDING EDGE AND GUIDE POINTS',I10/
     . '          NUMBER OF EDGE AND GUIDE POINTS ADDED   ',I10)
106     FORMAT(' X,Y,DX,DY,W   ',5F10.2)
107     FORMAT(' NUMBER OF CORRECTION VECTORS TO CALCULATE',2I10)
108     FORMAT(' Min number data points in ISTEPxISTEP box - below',
     . ' which a guide point will be added',I10)
150     FORMAT(/' SUM OF SQUARES OF RESIDUAL AT DATA POINTS=',G20.2)
151     FORMAT(' RANK OF SYSTEM    =',I10/
     . ' VALUE OF NCREAL USED  =',I10/
     . ' VALUE OF EPS USED =',F15.3)
159     FORMAT(': IMAXCOR ON INPUT =',I6,
     . '  , REDUCED TO BE EQUAL TO ISTEP =',I6)
160     FORMAT('$ ITYPE,IOUT,IMAXCOR,ISTEP ?')
161     FORMAT('               ITYPE----------',I5/
     . '               IOUT-----------',I5/
     . '               IMAXCOR--------',I5/
     . '               ISTEP----------',I5/
     . '               LTAPER---------',4X,L1/
     . '               RTAPER---------',F7.1/
     . '               LTABOUT--------',4X,L1)
162     FORMAT('$ FULL FILENAME OF OUTPUT IMAGE FILE ?')
163     FORMAT(A)
164     FORMAT(' IMAGE FILE CREATED WAS  ',A)
165     FORMAT('$ TITLE TO ADD TO OUTPUT IMAGE FILE ?')
166     FORMAT(20A4)
167     FORMAT(' TITLE ADDED TO IMAGE FILE WAS  ',20A4)
168     FORMAT('               IKNOTX---------',I5,/,
     . '               IKNOTY---------',I5,/,
     . '               EPS------------',F10.6,/,
     . '               FACTOR---------',F10.3,/,
     . '               THRESH---------',F10.1,/,
     . '               TLTAXIS--------',F10.1)
169     FORMAT('$ IKNOTX,IKNOTY,EPS,FACTOR,TLTAXIS ?')
170     FORMAT('::',I5,' PROGRAM DIMENSIONS INADEQUATE, PLEASE CHECK',/,
     . ' 1=ISIZEX;2=ISTEP;3=NDIMX;4=NDIMY;5=NPOINT(NDATA);',
     . ' 6=NCMAX;,7=NWSPCE;8=NPOINT(MCALC)',/)
171     FORMAT(':: PROGRAM DIMENSION IDEEP INADEQUATE, PLEASE CHECK') 
172     FORMAT(':: PRESENT PROG DIM INADEQUATE: MAX KNOTS ALLOWED=',I5)
191     FORMAT(':: ERROR *** IXA OUTSIDE RANGE 1 - ISIZEX',3I14,F15.6)
192     FORMAT(':: ERROR *** IYA OUTSIDE RANGE 1 - ITEST',4I14,F15.6)
193     FORMAT(' I,J,IX,IY,IXTRUE,IYTRUE=',12I6)
194     FORMAT(' T,XDELTA,YDELTA,XAPPLY,YAPPLY=',5F15.4)
195     FORMAT(' INTERPOLATION DONE, THIS STRIP NOW OUTPUT, J=',I6,
     . '  INOUT=',I5)
1101    FORMAT(':: FAILED KNOT SORTING IN OUTPUT E02ZAF, IFAIL=',I5)
1102    FORMAT(':: CUBIC SPLINE CALCULATION OF X-ERROR, IFAIL=',I5)
1103    FORMAT(':: CUBIC SPLINE CALCULATION OF Y-ERROR, IFAIL=',I5)
1105    FORMAT(' I,XPOS(I),YPOS(I),K,XOUT(K),YOUT(K)=',2(I10,2F10.2))
1106    FORMAT('  BELOW ARE ANY CORR PTS MORE THAN',F7.1,' FROM DATA')
3210    FORMAT(' TITLE FROM CCORDATA FILE; ',20A4)
3211    FORMAT(' FULL FILENAME OF INPUT IMAGE FILE? ')
3212    FORMAT('  FILENAME OF IMAGE FILE TO BE USED ',A)
8000    FORMAT(' TIME FOR STEP UP TO START OF BICUBIC SPLINE',F10.1)
8001    FORMAT(' TIME FOR BICUBIC SPLINE FITTING OF X-DIR',F10.1)
8002    FORMAT(' TIME FOR BICUBIC SPLINE FITTING OF Y-DIR',F10.1)
8003    FORMAT(' TIME FOR UNBENDING IMAGE',F10.1)
8004    FORMAT(' TIME UP TO START OF PLOT WITH ITYPE=0',2F10.1)
        END
C
C*******************************************************************************
C
C  SUBROUTINE TO PLOT THE LATTICE OF FITTED DISTORTION CORRECTIONS.
        SUBROUTINE PLOTCORR(TITLE,MCALC,XOUT,YOUT,XCORR,YCORR,NXYZ,IP,
     .   RMAG,LCOLOR)
        LOGICAL LCOLOR
        DIMENSION XOUT(1),YOUT(1),XCORR(1),YCORR(1)
        DIMENSION TITLE(20),NXYZ(3),TITLEPLOT(20),TEXT(20)
        CHARACTER DAT*24
CTSH++
        CHARACTER*80 TMPTITLEPLOT
        EQUIVALENCE (TMPTITLEPLOT,TITLEPLOT)
CTSH--
CHEN>
        DIMENSION line(20)
        CHARACTER*80 TMPline
        EQUIVALENCE (TMPline,line)
C
C-------Magnification (Exaggeration) factor 10 times:
C        RMAG = 10.0
CHEN<

        CALL FDATE(DAT)
        WRITE(6,1502) DAT(5:24)
1502    FORMAT('  Date from fdate ----  ',A20)
CTSH    WRITE(TITLEPLOT) (TITLE(J),J=1,15),DAT(5:24,1501)
CTSH++
        WRITE(TMPTITLEPLOT,1501) (TITLE(J),J=1,15),DAT(5:24)
CTSH--
1501    FORMAT(15A4,A20)
CHEN>
        WRITE(TMPline,'('' Vectors are plotted '',F4.1,'' times '',
     .    ''elongated.'')') RMAG
CHEN<
200     FORMAT('  ENTERING PLOTCORR')
201     FORMAT(20A4)
       WRITE(6,200)
       WRITE(6,201)TITLEPLOT
      PLTSIZ=260.0
      FONTSIZE=3.6
      IF(IP.EQ.100)CALL P2K_OUTFILE('CCPLOT.PS',9)
      CALL P2K_HOME
      CALL P2K_FONT('Courier'//CHAR(0),FONTSIZE)
      CALL P2K_GRID(0.5*PLTSIZ,0.5*PLTSIZ,1.0)
      CALL P2K_ORIGIN(-0.5*PLTSIZ,-0.7*PLTSIZ,0.)
      CALL P2K_COLOUR(0)
        SPLOT=PLTSIZ/MAX(NXYZ(1),NXYZ(2))
C  BOX ROUND THE WHOLE IMAGE AREA
                SIZEX=NXYZ(1)*SPLOT
                SIZEY=NXYZ(2)*SPLOT
      CALL P2K_MOVE(0.,0.,0.)
      CALL P2K_DRAW(SIZEX,0.,0.)
      CALL P2K_DRAW(SIZEX,SIZEY,0.)
      CALL P2K_DRAW(0.,SIZEY,0.)
      CALL P2K_DRAW(0.,0.,0.)
C NOW PLOT EACH ERROR VECTOR STARTING WITH A CROSS AT PROPER POSITION.
C  CROSS TEMPORARILY DISABLED 16.8.84
      CALL P2K_FONT('Courier'//CHAR(0),FONTSIZE*0.5)
CHEN>
      rlenmax=0.0
      DO 99  I=1,MCALC
        XPLOT=SPLOT*XOUT(I)
        YPLOT=SPLOT*YOUT(I)
        IF(XPLOT.GT.SIZEX.OR.XPLOT.LT.0.) GO TO 99
        IF(YPLOT.GT.SIZEY.OR.YPLOT.LT.0.) GO TO 99
        XERR=SPLOT*XCORR(I)*RMAG + XPLOT 
        YERR=SPLOT*YCORR(I)*RMAG + YPLOT 
        IF(XERR.GT.SIZEX.OR.XERR.LT.0.) GO TO 99
        IF(YERR.GT.SIZEY.OR.YERR.LT.0.) GO TO 99
        rlen=sqrt((XERR-XPLOT)**2+(YERR-YPLOT)**2)
        if(rlen.gt.rlenmax)rlenmax=rlen
99    CONTINUE
C-----Adjust maximal length to width/20
      rlenmax=rlenmax/20.0
      if(rlenmax.lt.0.001)rlenmax=0.001
CHEN<
      DO 100 I=1,MCALC
        XPLOT=SPLOT*XOUT(I)
        YPLOT=SPLOT*YOUT(I)
        IF(XPLOT.GT.SIZEX.OR.XPLOT.LT.0.) GO TO 100
        IF(YPLOT.GT.SIZEY.OR.YPLOT.LT.0.) GO TO 100
C       WRITE(TEXT,102)
C       CALL P2K_CSTRING(TEXT,1,0.)
102     FORMAT('X')
        XERR=SPLOT*XCORR(I)*RMAG + XPLOT
        YERR=SPLOT*YCORR(I)*RMAG + YPLOT
        IF(XERR.GT.SIZEX.OR.XERR.LT.0.) GO TO 100
        IF(YERR.GT.SIZEY.OR.YERR.LT.0.) GO TO 100
        CALL P2K_MOVE(XPLOT,YPLOT,0.)
CHEN>
C---------Calculate angle and length of line:
          rang=atan2(XERR-XPLOT,YERR-YPLOT)*180.0/3.141592654
          if(rang.lt.0.0)rang=rang+360.0
          rlen=sqrt((XERR-XPLOT)**2+(YERR-YPLOT)**2)/rlenmax
          rlen=(rlen*0.5)+0.5
          if(rlen.gt.1.0)rlen=1.0
          if(rlen.lt.0.5)rlen=0.5
C---------Assign HSV Values, as Hue=angle, Saturation=length, Value=1.0
          rhue=rang
          rsat=rlen
          rval=0.7
C---------Calculate RGB Values:
          c=rval * rsat
          rhue60=rhue / 60.0
          x = c * (1 - abs(mod(rhue60,2.0) - 1))
          if(rhue60<1.0)then
            rgbr=c
            rgbg=x
            rgbb=0
          elseif (rhue60<2.0)then
            rgbr=x
            rgbg=c
            rgbb=0
          elseif (rhue60<3.0)then
            rgbr=0
            rgbg=c
            rgbb=x
          elseif (rhue60<4.0)then
            rgbr=0
            rgbg=x
            rgbb=c
          elseif (rhue60<5.0)then
            rgbr=x
            rgbg=0
            rgbb=c
          else
            rgbr=c
            rgbg=0
            rgbb=x
          endif
          rgbr=rgbr+rval-c
          rgbg=rgbg+rval-c
          rgbb=rgbb+rval-c
C          write(*,'(''c,x='',2F6.3,'' rlen,rang='',2F9.3,'' r,g,b='',3F6.3)')
C     1      c,x,rlen,rang,rgbr,rgbg,rgbb
          if(rgbr.gt.1.0)rgbr=1.0
          if(rgbg.gt.1.0)rgbg=1.0
          if(rgbb.gt.1.0)rgbb=1.0
          if(rgbr.lt.0.0)rgbr=0.0
          if(rgbg.lt.0.0)rgbg=0.0
          if(rgbb.lt.0.0)rgbb=0.0
C---------Set RGB color:
          if(LCOLOR) CALL P2K_RGB_COLOUR(rgbr,rgbg,rgbb)
        CALL P2K_DRAW(XERR,YERR,0.)
CHEN<
100   CONTINUE
      CALL P2K_COLOUR(0)
      CALL P2K_FONT('Courier'//CHAR(0),FONTSIZE)
      CALL P2K_MOVE(10.,SIZEY+15.0,0.)
      CALL P2K_STRING(TITLEPLOT,80,0.)
      CALL P2K_MOVE(10.,SIZEY+5.0,0.)
      CALL P2K_STRING(line,80,0.)
                CALL P2K_PAGE
        RETURN
        END
C
C*******************************************************************************
C
      SUBROUTINE ROTATEXY(N,DX,DY,ANGLE)
      DIMENSION DX(1),DY(1)
      PI=3.1415926
      C=COS(ANGLE*PI/180.0)
      S=SIN(ANGLE*PI/180.0)
      DO 100 I = 1,N
        XNEW =  DX(I)*C + DY(I)*S
        YNEW = -DX(I)*S + DY(I)*C
      DX(I) = XNEW
100   DY(I) = YNEW
      RETURN
      END
C
C*******************************************************************************
C
      SUBROUTINE FILLEMPTIES(MXDIM,MYDIM,DXAV,DYAV,WAV,NAV)
CHENN>
C        PARAMETER (NDIMX=500)
C        PARAMETER (NDIMY=500)
        PARAMETER (NDIMX=2000)
        PARAMETER (NDIMY=2000)
CHENN<
        REAL DXAV(NDIMX,NDIMY),DYAV(NDIMX,NDIMY)        ! Bins for guide point
        REAL NAV(NDIMX,NDIMY),WAV(NDIMX,NDIMY)          ! and linear interpol.
                NZEROES=0
                DO 4000 I=1,MXDIM
                DO 4000 J=1,MYDIM
                        IF(NAV(I,J).EQ.0) NZEROES=NZEROES+1
4000            CONTINUE
                PROPZERO=FLOAT(NZEROES)/FLOAT(MXDIM*MYDIM)
                RATIO=SQRT(1.0/PROPZERO)
                IR=RATIO + 0.5
                NPUTIN=0
      DO 6000 ICYCLE=1,12
                ISTEP=2**(ICYCLE-1)
                IRC=IR*ISTEP
                DSTART=2.0*IRC**2
        DO 4300 I=1,MXDIM
        DO 4300 J=1,MYDIM
        IF(NAV(I,J).EQ.0) THEN  ! PUT IN NEAREST REAL CORRECTION
                DSQMIN=DSTART
                IDONE=0
                DO 4200 IS=I-IRC,I+IRC,ISTEP  ! inefficient but fairly simple
                DO 4200 JS=J-IRC,J+IRC,ISTEP
C                       IF(NAV(IS,JS).LE.0.OR.IS.LE.0.OR.JS.LE.0.OR.
C     .                         IS.GT.MXDIM.OR.JS.GT.MYDIM) GO TO 4200
CHENN>
C                       IF(IS.LE.0.OR.JS.LE.0.OR.IS.GT.MXDIM.OR.
C     *                         JS.GT.MYDIM.OR.NAV(IS,JS).LE.0.) GO TO 4200
                        IF(IS.LE.0.OR.JS.LE.0.OR.IS.GT.MXDIM) GO TO 4200
                        IF(JS.GT.MYDIM) GO TO 4200
                        IF(NAV(IS,JS).LE.0.) GO TO 4200
CHENN<
                        DSQ=(IS-I)**2+(JS-J)**2
                        IF(DSQ.LT.DSQMIN) THEN
                                IDONE=1
                                DSQMIN=DSQ
                                DXUSE=DXAV(IS,JS)
                                DYUSE=DYAV(IS,JS)
                                WUSE=WAV(IS,JS)/(1.0+DSQ)       ! w. low weight.
                        ENDIF
4200            CONTINUE
                IF(IDONE.EQ.1) THEN
                        NPUTIN=NPUTIN+1
                        DXAV(I,J)=DXUSE
                        DYAV(I,J)=DYUSE
                        WAV(I,J) =WUSE
                        NAV(I,J) = -1           ! distinguish new vectors from old
                ENDIF
        ENDIF
4300    CONTINUE
                DO 4310 I=1,MXDIM
                DO 4310 J=1,MYDIM
                        IF(NAV(I,J).LT.0) NAV(I,J)=-NAV(I,J)    ! remove distinction
4310            CONTINUE
                WRITE(6,4301) MXDIM*MYDIM,NZEROES,ICYCLE,
     . NPUTIN,IRC,NZEROES-NPUTIN
4301            FORMAT( ' Total number of bins         ',I7/
     . ' Number of empty bins         ',I7/
     . ' Number filled in cycle',I3,' was',
     . I7,' using range',I5/
     . ' Number still unfilled       ',I7/)
        IF(NZEROES-NPUTIN.EQ.0) RETURN
6000  CONTINUE
      RETURN
      END
C
C******************************************************************************
      SUBROUTINE CALCTAPER(ISTEP,TAPER,RTAPER)
C       apply a cosine bell taperedge function starting at radius RTAPER
C       and extending out to edge of ISTEP box.
      PARAMETER (ISTEPMAX=60)
      REAL TAPER(ISTEPMAX,ISTEPMAX)
      XYRAD = (ISTEP-1)/2.0
      IF(RTAPER.GT.XYRAD-1) RTAPER=XYRAD-1      ! at least one pixel clear space
      DO 65 I = 1,ISTEP
        DO 65 J = 1,ISTEP
          RAD = (I-1-XYRAD)**2 + (J-1-XYRAD)**2
          RAD = SQRT(RAD)
          IF(RAD.LE.RTAPER) THEN
           TAPER(I,J) = 1.0
          ELSE
            IF(RAD.GE.XYRAD-1) THEN
              TAPER(I,J) = 0.0
            ELSE
              ARG = (RAD-RTAPER)/(XYRAD-1-RTAPER)*0.5*3.1415926
              TAPER(I,J) = COS(ARG)
            ENDIF
          ENDIF
65    CONTINUE
      RETURN
      END
C
C*******************************************************************************
C
C  WRITE OUT THE DISTORTION CORRECTION TABLE AS USED IN THE UNBENDING 
C   PROGRAM SO THAT IT MIGHT BE USED IN ANOTHER PROGRAM.
C
        SUBROUTINE  WRITETABLE(TITLE,ISTEP,MXDIM,MYDIM,NXYZ,XCORR,YCORR)
        DIMENSION XCORR(1),YCORR(1)
        INTEGER NXYZ(3)
        CHARACTER*4 TITLE(40),TITLEOUT(40)
        INTEGER ISTEP,MXDIM,MYDIM
CHENN>
        INTEGER*8 IBIN(-1001:1001)
CHENN<
        CHARACTER DAT*24
CTSH++
        CHARACTER*80 TMPTITLEOUT
        EQUIVALENCE (TMPTITLEOUT,TITLEOUT)
CTSH--

                CALL CCPDPN(9,'TABLEOUT','UNKNOWN','F',0,0)
                CALL FDATE(DAT)
                WRITE(6,10) DAT(5:24)
10              FORMAT('  Date from fdate ----  ',A20)
CTSH                WRITE(TITLEOUT) (TITLE(J),J=1,15),DAT(5:24,11)
CTSH++
                WRITE(TMPTITLEOUT,11) (TITLE(J),J=1,15),DAT(5:24)
CTSH--
11              FORMAT(15A4,A20)
20      FORMAT('  ENTERING WRITETABLE')
21      FORMAT(20A4)
22      FORMAT(' ISTEP,MXDIM,MYDIM,NXYZ =',6I6)
      WRITE(6,20)
      WRITE(9,21) TITLEOUT
      WRITE(9,22) ISTEP,MXDIM,MYDIM,(NXYZ(J),J=1,3)
      MCALC=MXDIM*MYDIM
C Steps along X are indexed most rapidly in arrays XCORR,YCORR
CHENN>
C       DO 100 I=1,MCALC
C100            WRITE(9,30) XCORR(I),YCORR(I)
C30             FORMAT(2F8.2)
C
C        WRITE(6,*) XCORR(I),YCORR(I)
C        WRITE(9,*) XCORR(I),YCORR(I)
C          WRITE(6,30) XCORR(I)           ,YCORR(I),
C     1                XCORR(I)-XCORR(I-1),YCORR(I)-YCORR(I-1)
C           
C          WRITE(9,30) XCORR(I)           ,YCORR(I),
C     1                XCORR(I)-XCORR(I-1),YCORR(I)-YCORR(I-1)
100     continue
30      FORMAT(2F8.2,2F9.2)
CHENN<
        CLOSE(9)
CHENN>
C
        DO I=-1000,1000
          IBIN(I)=0.0
        enddo
        do I=2,MCALC
          xdiff=XCORR(I)-XCORR(I-1)
          ydiff=YCORR(I)-YCORR(I-1)
          idiff=INT(sqrt(xdiff*xdiff+ydiff*ydiff)*10)
          if(idiff.gt. 1000)idiff= 1000
          if(idiff.lt.-1000)idiff=-1000
          IBIN(idiff)=IBIN(idiff)+1
        enddo
        do i=-1000,1000
          rval=i/10.0
          write(17,'(F9.3,I15)')rval,IBIN(i)
        enddo
C
CHENN<
        RETURN
        END
C
C
c==========================================================
c
      SUBROUTINE shorten(czeile,k)
C
C counts the number of actual characters not ' ' in czeile
C and gives the result out in k.
C
      CHARACTER * (*) CZEILE
      CHARACTER * 1 CTMP1
      CHARACTER * 1 CTMP2
      CTMP2=' '
C
      ilen=len(czeile)
      DO 100 I=1,ilen
         k=ilen+1-I
         READ(CZEILE(k:k),'(A1)')CTMP1
         IF(CTMP1.NE.CTMP2)GOTO 300
  100 CONTINUE
  300 CONTINUE
      IF(k.LT.1)k=1
C
      RETURN
      END

